########################################################################
# R-code:	CASA_Stn_Vario.R
# Purpose:	For fitting variogram model paramteters with 'Run_vario.R',
# 			and kriging with 'Run_krige.r'.
#			Based upon anomalies of Central America/South America
#			rainfall station data.
########################################################################


############################
# LOAD REQUIRED 'R' PACKAGES

#help(available.packages)
#search()	#Lists currently loaded packages

#Might help with installing 'rgdal' or its dependency 'proj4'
#setRepositories(ind=1:2) 

#install.packages("sp", repos="http://cran.us.r-project.org")#
#install.packages("zoo", repos="http://cran.us.r-project.org")#
#install.packages("xts", repos="http://cran.us.r-project.org")#
#install.packages("spacetime", repos="http://cran.us.r-project.org")#
#install.packages("gstat", repos="http://cran.us.r-project.org")#
#install.packages("proj4", repos="http://cran.us.r-project.org")
#install.packages("rgdal", repos="http://cran.us.r-project.org")
#install.packages("raster", repos="http://cran.us.r-project.org")
#install.packages("rasterVis", repos="http://cran.us.r-project.org")
#install.packages("geoR", repos="http://cran.us.r-project.org")
#install.packages("MASS", repos="http://cran.us.r-project.org")
#install.packages("maps", repos="http://cran.us.r-project.org")
#install.packages("maptools", repos="http://cran.us.r-project.org")
#install.packages("mapdata", repos="http://cran.us.r-project.org")
#install.packages("splancs", repos="http://cran.us.r-project.org")
#install.packages("RandomFields", repos="http://cran.us.r-project.org")
#install.packages("geoR", repos="http://cran.us.r-project.org")
#install.packages("automap", repos="http://cran.us.r-project.org")
#install.packages("RColorBrewer", repos="http://cran.us.r-project.org")


library (sp)        	# sp        Depends: R (= 2.10.0), methods, graphics                   
library (zoo)       	# zoo       Depends: R (= 2.10.0), stats                      
library (xts)       	# xts       Depends: zoo (= 1.7-2)                    
library (spacetime) 	# spacetime Depends: R (= 2.10.0), methods, sp (= 0.9-89), zoo, xts                          
library (gstat)    	 	# gstat     Depends: R (= 2.10), methods, sp (= 0.9-72), spacetime (= 0.4-4) 
library (proj4)     	# proj4                        
library (rgdal)     	# rgdal     Depends: R (= 2.10.0), methods, sp (= 0.9-60)                        
library (raster)    	# raster    Depends: methods, sp, R (= 2.11.0) 
library (rasterVis)    	# rasterVis 
library (geoR)   		# geoR
library (MASS)   		# MASS
library (maps)  		# maps
library (maptools) 		# maptools
library (mapdata) 		# mapdata
library (splancs)		# splancs
library (RandomFields)	# RandomFields
library (geoR)			# geoR
library (automap)		# automap
library (RColorBrewer)  # RColorBrewer


#######
# PATHS
#
# Note that subdirectories require different slashes per operating system
# (i.e. "/" for Unix, "\\" for Windows)

DirPrecipSourceVgmKrg	=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\CODE\\"

DirRefData				=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\DATA\\FTIP\\"

DirLandMask				=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\DATA\\LandMask\\"


DirPrecipVgmPar			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\RESULTS\\Precip\\Variogram_Parameters\\"
DirPrecipData			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\DATA\\Stations\\"
DirPrecipMeta			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\DATA\\Stations\\"
DirPrecipKrige			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\RESULTS\\Precip\\"	# Note that subdirectories will be created for months
DirPrecipMosaic			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\RESULTS\\Precip\\"	# Note that subdirectories will be used for months
DirPrecipMosaicMask		=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\RESULTS\\Precip\\"	# Note that subdirectories will be used for months

#DirTempVgmPar			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\RESULTS\\Temp\\Variogram_Parameters\\"
#DirTempData			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\DATA\\Temp\\2012-05-02\\"
#DirTempMeta			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\DATA\\Temp\\2012-05-02\\"
#DirTempKrige			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\RESULTS\\Temp\\"	# Note that subdirectories will be created for months
#DirTempMosaic			=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\RESULTS\\Temp\\"	# Note that subdirectories will be used for months
#DirTempMosaicMask		=	"C:\\Users\\User\\git\\fews-02-casa\\fews-02-casa\\RESULTS\\Temp\\"	# Note that subdirectories will be used for months




FilePrecipSourceVgmKrg		=	"CASA_Stn_Vario.R"

FilePrecipVgmPar			=	"Variogram_Parameters.csv"	# Could use list.files for many files, but would need to update dependent code.
FileTempVgmPar				=	"Variogram_Parameters.csv"	# Could use list.files for many files, but would need to update dependent code.

FilePrefixPrecipRefData		=	"ftip2009"
FileSuffixPrecipRefData		=	".bil"

FileLandMask				=	"global_landmask.50N50S.p25Deg.tif"

FilePrefixPrecipData		=	"cacol"
FileSuffixPrecipData		=	"_mth.csv"
#FilePrefixPrecipMeta		=	"precp"
#FileSuffixPrecipMeta		=	"meta.csv"

FilePrefixPrecipKrgEstMm	=	"krmm"
FilePrefixPrecipKrgErrMm	=	"semm"
FilePrefixPrecipKrgEstPc	=	"krpc"
FilePrefixPrecipKrgErrPc	=	"sepc"
FileSuffixPrecipKrgImg		=	".bil"
FileSuffixPrecipKrgHdr		=	".hdr"

#FilePrefixTempData			=	"temp"
#FileSuffixTempData			=	"data.csv"
#FilePrefixTempMeta			=	"temp"
#FileSuffixTempMeta			=	"meta.csv"
#
#FilePrefixTempKrgEstDg		=	"krdeg"
#FilePrefixTempKrgErrDg		=	"sedeg"
#FileSuffixTempKrgImg		=	".bil"
#FileSuffixTempKrgHdr		=	".hdr"



##################
# GLOBAL VARIABLES

MinLatBG     	=  -50   	# Minimum latitude of background climatology (possibly specify with switch)
MaxLatBG     	=   50   	# Maximum latitude of background climatology (possibly specify with switch)
MinLonBG     	= -180  	# Minimum longitude of background climatology (possibly specify with switch)
MaxLonBG     	=  180  	# Maximum longitude of background climatology (possibly specify with switch)
MinValBG 		=    0		#
MaxValBG		= 5000		#

YearBegData  	= 2009  	# Beginning year of station data
YearEndData  	= 2009  	# End year of station data
YearBegRefMn  	= 2009  	# Beginning year of refernce mean for anomaly calculation
YearEndRefMn   	= 2009  	# End year of refernce mean for anomaly calculation
YearBegVgm  	= 2009  	# Beginning year of variogram
YearEndVgm   	= 2009  	# End year of variogram
YearBegKrg  	= 2009  	# Beginning year of kriging
YearEndKrg   	= 2009  	# End year of kriging

VgmModelForm 	= "Exp"		# Covariogram model form (as text, i.e. "Exp", "Sph", etc.)
MaxVarioDist 	=   14		# Maximum distance for calculating variogram model

BufferReg		=   2.		# Buffer around region to address edge effects by including nearby stations.

SigmaPrecip		=    4		# Sigma probability value for removing extreme values from precipitation data.
SigmaTemp		=    4		# Sigma probability value for removing extreme values from temperature data.

KrigeNMax		=   20		# Maximum number of nearest observations used for Kriging

CellLatDim		= 0.05		# Dimension length of latitude for each raster cell
CellLonDim		= 0.05		# Dimension length of longitude for each raster cell

WinLatDimPre	=    1		# Dimension length of latitude before  focal cell for moving window variogram parameter estimation
WinLatDimPost	=    1		# Dimension length of latitude after   focal cell for moving window variogram parameter estimation
WinLonDimPre	=    1		# Dimension length of longitude before focal cell for moving window variogram parameter estimation
WinLonDimPost	=    1		# Dimension length of longitude after  focal cell for moving window variogram parameter estimation

#LatRegionDim	=   50		# Dimension length of latitude for subdividing entire region of interest into smaller regions
#LatMinNominal	=  -50		# Minimum latitude value of entire region of interest (as named boundary edge, later processed for raster cell center)
#LatMaxNominal	=   50		# Maximum latitude value of entire region of interest (as named boundary edge, later processed for raster cell center)
#LonRegionDim	=   60		# Dimension length of longitude for subdividing entire region of interest into smaller regions
#LonMinNominal	= -180		# Minimum longitude value of entire region of interest (as named boundary edge, later processed for raster cell center)
#LonMaxNominal	=  180		# Maximum longitude value of entire region of interest (as named boundary edge, later processed for raster cell center)




###########
# FILE INFO

# Header column names of input data
#DataFileColNames	= seq.int(YearBegData, YearEndData)
DataFileColNames	= c("Sts_ID","Longitutde","Latitude","year","mth1","mth2","mth3","mth4","mth5","mth6","mth7","mth8","mth9","mth10","mth11","mth12") 


## Header column names of input metadata
#MetaFileColNames	= c("CountryCode","StnName","Lat","Lon","Elev","StnSeqNum","NativeID")




##################
# PRELIMINARY CODE

# New for CASA

#source(paste(DirPrecipSourceVgmKrg, FilePrecipSourceVgmKrg, sep=""))

 year = 09
#for (year in 09:09){
	
	#Read data (CSV file per year with columns for station ID, year, long, lat, and monthly rainfall)
	dataStn = 	read.csv(paste(DirPrecipData, FilePrefixPrecipData, sprintf("%02d", year), FileSuffixPrecipData,  sep=""),
						header=TRUE, stringsAsFactors=FALSE, strip.white=TRUE)	
				
	# Rounding values for region, though could use original values
	LatMinNominal	=	floor	(min(dataStn$Latitude))			# Minimum latitude value of entire region of interest (as named boundary edge, later processed for raster cell center)
	LatMaxNominal	=	ceiling	(max(dataStn$Latitude))			# Maximum latitude value of entire region of interest (as named boundary edge, later processed for raster cell center)
	LatRegionDim	=	abs		(LatMaxNominal-LatMinNominal)	# Dimension length of latitude for subdividing entire region of interest into smaller regions
	LonMinNominal	=	floor	(min(dataStn$Longitude))		# Minimum longitude value of entire region of interest (as named boundary edge, later processed for raster cell center)
	LonMaxNominal	=	ceiling	(max(dataStn$Longitude))		# Maximum longitude value of entire region of interest (as named boundary edge, later processed for raster cell center)
	LonRegionDim	=	abs		(LonMaxNominal-LonMinNominal)	# Dimension length of longitude for subdividing entire region of interest into smaller regions
	
#	defineRegions()
	
	
	 month = 8
#	for (month in 1:12)	{
		
		outputDir = sprintf("%02d", month)
		dir.create(file.path(DirPrecipVgmPar, outputDir), recursive=TRUE, showWarnings=FALSE)
		outputDir = paste (outputDir, "\\", sep="")
		colName = paste("mth", month, sep="")
		colNumber = which(colnames(dataStn)==colName)
		dataStnMon	= data.frame(dataStn$Longitude, dataStn$Latitude, dataStn[,colNumber])
		colnames(dataStnMon) <- c("Longitude","Latitude","Rain_Stn")
		
		
		# Find indices of data that have
		# values below 0 and above 900
		# for assignment of NA values.
		ok  = dataStnMon$Rain_Stn > 0 # | dataStnMon < 900
		
		# Matrix with data and NA values as specified
#		dataStnMon[setNA]   = NA
		dataStnMon	= data.frame(dataStnMon[ok,])
		
		dedup = function(data) { 							# "deduplicate": remove duplicate lon, lat values
			
			lon  = data[,1]
			lat  = data[,2]
			ll   = cbind(lon,lat)
			lln = unique(ll)								# return only unique lon, lat data
			dn  = matrix(NA,nrow=nrow(lln),ncol=ncol(data))	# prepare matrix for output
			colnames(dn) = colnames(data)
			nc = ncol(data)
			for( i in 1:nrow(dn) ) {
				sub = lon == lln[i,1] & lat == lln[i,2]		#identify original rows (stations) with the same lon, lat values as the unique lon, lat
				for( j in 1:nc )
					dn[i,j] = mean(data[sub,j])				#for each column, place the mean of each row or duplicate sets of rows into ouput data
			}
			return(dn)
		}		
		
		# De-duplicate (find and average duplicates)
		dataStnMon	= data.frame(dedup(dataStnMon))
		

		dataStnMonExt		= extent( c(xmin= floor		(min	(dataStnMon$Longitude)) - BufferReg, 
										xmax= ceiling	(max	(dataStnMon$Longitude)) + BufferReg, 
										ymin= floor		(min	(dataStnMon$Latitude))  - BufferReg, 
										ymax= ceiling	(max	(dataStnMon$Latitude))  + BufferReg))
										
		smallTestExt		= extent( c(xmin= -80, 
										xmax= -75, 
										ymin= 5, 
										ymax= 10))
				
				
		FTIPMonImg			= raster(paste(DirRefData, FilePrefixPrecipRefData, sprintf("%02d", month), FileSuffixPrecipRefData, sep=""), 
									crs= "+proj=lonlat +datum=WGS84")
						
		FTIPMonImg			= crop(FTIPMonImg, dataStnMonExt)
		
		FTIPMonSmallTest	= crop(FTIPMonImg, smallTestExt)
		
		Rain_Sat			= extract(FTIPMonImg, as.data.frame(cbind(dataStnMon$Longitude,dataStnMon$Latitude)))

	
		
		Rain_Diff			= Rain_Sat - dataStnMon$Rain_Stn
	
		dataStnMon 			= as.data.frame(cbind(dataStnMon, Rain_Sat, Rain_Diff))

		
		
		## Stats Stn
		dataStnMonMax 		= max		(dataStnMon$Rain_Stn, na.rm = TRUE)
		dataStnMonMean 		= mean		(dataStnMon$Rain_Stn, na.rm = TRUE)
		dataStnMonMedian 	= median	(dataStnMon$Rain_Stn, na.rm = TRUE)
		dataStnMonMin 		= min		(dataStnMon$Rain_Stn, na.rm = TRUE)
		#dataStnMonVar 		= var		(dataStnMon$Rain_Stn, na.rm = TRUE)
		dataStnMonSD 		= sd		(dataStnMon$Rain_Stn, na.rm = TRUE)
		
		cat(paste(	#"Region: ... (if looped)
						"Year:      ", round(year), 			"\n",
						"Month:     ", round(month), 			"\n",
						"Max:       ", round(dataStnMonMax), 	"\n",
						"Mean:      ", round(dataStnMonMean),	"\n", 
						"Median:    ", round(dataStnMonMedian),	"\n",
						"Min:       ", round(dataStnMonMin), 	"\n",
						#"Variance:  ", round(dataStnMonVar), 	"\n",
						"Std. Dev.: ", round(dataStnMonSD), 	"\n",
						sep = ""),
		file = paste(DirPrecipVgmPar, outputDir, "BasicStats_Stn_", sprintf("%02d", month), ".txt", sep=""))
				
		## Stats Sat
		dataSatMonMax 		= max		(dataStnMon$Rain_Sat, na.rm = TRUE)
		dataSatMonMean 		= mean		(dataStnMon$Rain_Sat, na.rm = TRUE)
		dataSatMonMedian 	= median	(dataStnMon$Rain_Sat, na.rm = TRUE)
		dataSatMonMin 		= min		(dataStnMon$Rain_Sat, na.rm = TRUE)
		#dataSatMonVar 		= var		(dataStnMon$Rain_Sat, na.rm = TRUE)
		dataSatMonSD 		= sd		(dataStnMon$Rain_Sat, na.rm = TRUE)
		
		cat(paste(	#"Region: ... (if looped)
						"Year:      ", round(year), 			"\n",
						"Month:     ", round(month), 			"\n",
						"Max:       ", round(dataSatMonMax), 	"\n",
						"Mean:      ", round(dataSatMonMean),	"\n", 
						"Median:    ", round(dataSatMonMedian),	"\n",
						"Min:       ", round(dataSatMonMin), 	"\n",
						#"Variance:  ", round(dataSatMonVar), 	"\n",
						"Std. Dev.: ", round(dataSatMonSD), 	"\n",
						sep = ""),
				file = paste(DirPrecipVgmPar, outputDir, "BasicStats_Sat_", sprintf("%02d", month), ".txt", sep=""))

		## Stats Diff
		dataDiffMonMax 		= max		(dataStnMon$Rain_Diff, na.rm = TRUE)
		dataDiffMonMean 	= mean		(dataStnMon$Rain_Diff, na.rm = TRUE)
		dataDiffMonMedian 	= median	(dataStnMon$Rain_Diff, na.rm = TRUE)
		dataDiffMonMin 		= min		(dataStnMon$Rain_Diff, na.rm = TRUE)
		#dataDiffMonVar 	= var		(dataStnMon$Rain_Diff, na.rm = TRUE)
		dataDiffMonSD 		= sd		(dataStnMon$Rain_Diff, na.rm = TRUE)
		
		cat(paste(	#"Region: ... (if looped)
						"Year:      ", round(year), 			"\n",
						"Month:     ", round(month), 			"\n",
						"Max:       ", round(dataDiffMonMax), 	"\n",
						"Mean:      ", round(dataDiffMonMean),	"\n", 
						"Median:    ", round(dataDiffMonMedian),	"\n",
						"Min:       ", round(dataDiffMonMin), 	"\n",
						#"Variance:  ", round(dataDiffMonVar), 	"\n",
						"Std. Dev.: ", round(dataDiffMonSD), 	"\n",
						sep = ""),
				file = paste(DirPrecipVgmPar, outputDir, "BasicStats_Diff_", sprintf("%02d", month), ".txt", sep=""))

		
		
		hist(dataStnMon$Rain_Stn, main=paste("Histogram of (Station) Rainfall Month ", sprintf("%02d", month), sep=""), xlab="Rainfall (mm)")
		dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Histogram_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
		
		hist(dataStnMon$Rain_Sat, main=paste("Histogram of (Satellite: FTIP) Rainfall Month ", sprintf("%02d", month), sep=""), xlab="Rainfall (mm)")
		dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Histogram_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
		
		hist(dataStnMon$Rain_Diff, main=paste("Histogram of (Station-Satellite Diff.) Rainfall Month ", sprintf("%02d", month), sep=""), xlab="Rainfall (mm)")
		dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Histogram_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
		
		



		# Plot Rain_Stn
		plot	(dataStnMon$Longitude, dataStnMon$Latitude, 
				col=colorRampPalette(c("brown","gray80","green"))(5)[cut(log(dataStnMon$Rain_Stn),5)],
				pch=19, cex=0.75, 
				main= "Rainfall: Station Observations",
				xlab= "Longitude", ylab= "Latitude")
		map		('worldHires', add=TRUE, fill= TRUE, col="gray98",
				xlim= range(dataStnMon$Longitude), 
				ylim= range(dataStnMon$Latitude))
		points	(dataStnMon$Longitude, dataStnMon$Latitude, 
				col=colorRampPalette(c("brown","gray80","green"))(5)[cut(log(dataStnMon$Rain_Stn),5)],
				pch=19, cex=0.75)
		legend	("bottomleft", 
				legend=c(paste(round(max(dataStnMon$Rain_Stn)), " mm", sep=""), "0 mm"), 
				col=c("green", "brown"), pch=19, cex=0.75)		
		dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Rain_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
		
		
		# Plot Rain_Sat
		plot	(dataStnMon$Longitude, dataStnMon$Latitude, 
				col=colorRampPalette(c("brown","gray80","green"))(5)[cut(log(dataStnMon$Rain_Sat),5)],
				pch=19, cex=0.75, 
				main= "Rainfall: Satellite (FTIP) Observations",
				xlab= "Longitude", ylab= "Latitude")
		map		('worldHires', add=TRUE, fill= TRUE, col="gray98",
				xlim= range(dataStnMon$Longitude), 
				ylim= range(dataStnMon$Latitude))
		points	(dataStnMon$Longitude, dataStnMon$Latitude, 
				col=colorRampPalette(c("brown","gray80","green"))(5)[cut(log(dataStnMon$Rain_Sat),5)],
				pch=19, cex=0.75)
		legend	("bottomleft", 
				legend=c(paste(round(max(dataStnMon$Rain_Sat)), " mm", sep=""), "0 mm"), 
				col=c("green", "brown"), pch=19, cex=0.75)				
		dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

		
		# Plot Rain_Diff
		plot	(dataStnMon$Longitude, dataStnMon$Latitude, 
				col=colorRampPalette(c("brown","gray80","green"))(5)[cut(dataStnMon$Rain_Diff,5)],
				pch=19, cex=0.75, 
				main= "Rainfall: (Station - Satellite) Differences",
				xlab= "Longitude", ylab= "Latitude")
		map		('worldHires', add=TRUE, fill= TRUE, col="gray98",
				xlim= range(dataStnMon$Longitude), 
				ylim= range(dataStnMon$Latitude))
		points	(dataStnMon$Longitude, dataStnMon$Latitude, 
				col=colorRampPalette(c("brown","gray80","green"))(5)[cut(dataStnMon$Rain_Diff,5)],
				pch=19, cex=0.75)
		legend	("bottomleft", 
				legend=c(paste(round(max(dataStnMon$Rain_Diff)), " mm", sep=""), paste(round(min(dataStnMon$Rain_Diff)), " mm", sep="")), 
				col=c("green", "brown"), pch=19, cex=0.75)		
		dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Rain_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

		
												

#########################
#		geoR package:
#
#		variog							calculates sample variogram
#		variog4 						directional variograms
#		variog.model.env 				variogram envelopes
#		variog.mc.env					variogram envelopes
#		variofit						variogram based parameter estimation 
#		plot.variogram					graphical output 

dataStnMonGeo 			= as.geodata(dataStnMon, coords.col = 1:2, data.col = 3:5, data.names = c("Rain_Stn", "Rain_Sat", "Rain_Diff"))

dataStnMonGeoVgmStn		= variog(dataStnMonGeo, coords = dataStnMonGeo$coords, data = dataStnMonGeo$data[,1], max.dist=3)
dataStnMonGeoVgmSat		= variog(dataStnMonGeo, coords = dataStnMonGeo$coords, data = dataStnMonGeo$data[,2], max.dist=3)
dataStnMonGeoVgmDiff	= variog(dataStnMonGeo, coords = dataStnMonGeo$coords, data = dataStnMonGeo$data[,3], max.dist=3)

plot(dataStnMonGeoVgmStn, main="[geoR variog] Variogram: Stations")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_geoR_Rain_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
plot(dataStnMonGeoVgmSat, main="[geoR variog] Variogram: Satellite (FTIP)")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_geoR_Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
plot(dataStnMonGeoVgmDiff, main="[geoR variog] Variogram: (Station - Satellite) Differences")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_geoR_Rain_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)



dataStnMonGeoVgmModFitStn		= variofit(dataStnMonGeoVgmStn, ini.cov.pars=c(27000, 1))
dataStnMonGeoVgmModFitSat		= variofit(dataStnMonGeoVgmSat, ini.cov.pars=c(12000, 1))
dataStnMonGeoVgmModFitDiff		= variofit(dataStnMonGeoVgmDiff, ini.cov.pars=c(11000, 1))

plot(dataStnMonGeoVgmStn, main="[geoR variofit] Variogram: Stations")
lines(dataStnMonGeoVgmModFitStn)
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_geoR_fit_Rain_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
plot(dataStnMonGeoVgmSat, main="[geoR variofit] Variogram: Satellite (FTIP)")
lines(dataStnMonGeoVgmModFitSat)
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_geoR_fit_Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
plot(dataStnMonGeoVgmDiff, main="[geoR variofit] Variogram: (Station - Satellite) Differences")
lines(dataStnMonGeoVgmModFitDiff)
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_geoR_fit_Rain_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)






#########################
#		'gstat' package:
#
#		variogram						Calculates the sample variogram from data
#		vgm								Generates a variogram model, or adds to an existing model
#		fit.variogram					Fit ranges and/or sills from a simple or nested variogram model

dataStnMonStn			= dataStnMon[,c(1,2,3)]
dataStnMonStnGstat 		<- gstat(formula = Rain_Stn~1, locations = ~Longitude+Latitude, data = dataStnMonStn)

dataStnMonSat			= dataStnMon[,c(1,2,4)]
dataStnMonSatGstat 		<- gstat(formula = Rain_Sat~1, locations = ~Longitude+Latitude, data = dataStnMonSat)

dataStnMonDiff			= dataStnMon[,c(1,2,5)]
dataStnMonDiffGstat 	<- gstat(formula = Rain_Diff~1, locations = ~Longitude+Latitude, data = dataStnMonDiff)


dataStnMonStnGstatVgm		= variogram(dataStnMonStnGstat, cutoff=3)
plot(dataStnMonStnGstatVgm, col="black", main="[gstat variogram] Variogram: Stations")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_gstat_Rain_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

dataStnMonSatGstatVgm		= variogram(dataStnMonSatGstat, cutoff=3)
plot(dataStnMonSatGstatVgm, col="black", main="[gstat variogram] Variogram: Satellite (FTIP)")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_gstat_Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

dataStnMonDiffGstatVgm		= variogram(dataStnMonDiffGstat, cutoff=3)
plot(dataStnMonDiffGstatVgm, col="black", main="[gstat variogram] Variogram: (Station - Satellite) Differences")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_gstat_Rain_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)




dataStnMonStnGstatVgmModFit	= fit.variogram(dataStnMonStnGstatVgm,  vgm(27000,"Exp",1,0))
plot(dataStnMonStnGstatVgm, model=dataStnMonStnGstatVgmModFit, col="black", main="[gstat fit.variogram] Variogram: Stations")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_gstat_fit_Rain_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

dataStnMonSatGstatVgmModFit	= fit.variogram(dataStnMonSatGstatVgm,  vgm(12000,"Exp",1,0))
plot(dataStnMonSatGstatVgm, model=dataStnMonSatGstatVgmModFit, col="black", main="[gstat fit.variogram] Variogram: Satellite (FTIP)")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_gstat_fit_Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

dataStnMonDiffGstatVgmModFit	= fit.variogram(dataStnMonDiffGstatVgm,  vgm(11000,"Exp",1,0))
plot(dataStnMonDiffGstatVgm, model=dataStnMonDiffGstatVgmModFit, col="black", main="[gstat fit.variogram] Variogram: (Station - Satellite) Differences")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_gstat_fit_Rain_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)






#########################
#		'automap' package:
#
#		autofitVariogram				Automatically fitting a variogram
#		autoKrige						Performs an automatic interpolation
#		autoKrige.cv					Automatic cross-validation
#		automapPlot						Special plot function for automap
#		compare.cv						Comparing the results of cross-validations
#		plot.autoKrige					Plot methods in automap
#		plot.posPredictionInterval		Plot methods in automap
#		posPredictionInterval			Determines the position of the p% prediction interval
		
#StnMonSpDataVgmAutoFit	= autofitVariogram	(StnMonSpData@data~1, locations= ~ StnMonSpData@coords[,1] + StnMonSpData@coords[,2], StnMonSpData@data, cutoff=10, width= 1, map=TRUE)
#plot(StnMonSpDataVgmAutoFit)

dataStnMonSpdf 					= SpatialPointsDataFrame(coords= dataStnMon[,1:2], data = dataStnMon[,3:5])

dataStnMonStnSpdf 				= SpatialPointsDataFrame(coords= dataStnMon[,1:2], data = as.data.frame(dataStnMon$Rain_Stn))
names(dataStnMonStnSpdf@data) 	<- "Rain_Stn"

dataStnMonSatSpdf 				= SpatialPointsDataFrame(coords= dataStnMon[,1:2], data = as.data.frame(dataStnMon$Rain_Sat))
names(dataStnMonSatSpdf@data) 	<- "Rain_Sat"

dataStnMonDiffSpdf 				= SpatialPointsDataFrame(coords= dataStnMon[,1:2], data = as.data.frame(dataStnMon$Rain_Diff))
names(dataStnMonDiffSpdf@data) 	<- "Rain_Diff"




dataStnMonStnAutoVgmModFit		= autofitVariogram(formula= 	dataStnMonStnSpdf@data[,1] ~ 1, 
													input_data= dataStnMonStnSpdf,
													verbose=	TRUE)										
plot(dataStnMonStnGstatVgm, 
		model=dataStnMonStnAutoVgmModFit$var_model, 
		main="[automap autofitVariogram] Variogram: Stations",
		sub=paste("Model:", 		dataStnMonStnAutoVgmModFit$var_model[2,1], 
				",   Nugget: ", 	round(dataStnMonStnAutoVgmModFit$var_model[1,2]), 
				",   Full Sill: ", 	round(dataStnMonStnAutoVgmModFit$var_model[2,2] + dataStnMonStnAutoVgmModFit$var_model[1,2]), 
				",   Range: ", 		round(dataStnMonStnAutoVgmModFit$var_model[2,3], digits=2), 
				sep=""),
		col="black")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_automap_fit_Rain_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)



dataStnMonSatAutoVgmModFit		= autofitVariogram(formula= 	dataStnMonSatSpdf@data[,1] ~ 1, 
		input_data= dataStnMonSatSpdf,
		verbose=	TRUE)
plot(dataStnMonSatGstatVgm, 
		model=dataStnMonSatAutoVgmModFit$var_model, 
		main="[automap autofitVariogram] Variogram: Satellite (FTIP)",
		sub=paste("Model:", 		dataStnMonSatAutoVgmModFit$var_model[2,1], 
				",   Nugget: ", 	round(dataStnMonSatAutoVgmModFit$var_model[1,2]), 
				",   Full Sill: ", 	round(dataStnMonSatAutoVgmModFit$var_model[2,2] + dataStnMonSatAutoVgmModFit$var_model[1,2]), 
				",   Range: ", 		round(dataStnMonSatAutoVgmModFit$var_model[2,3], digits=2), 
				sep=""),
		col="black")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_automap_fit_Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)



dataStnMonDiffAutoVgmModFit		= autofitVariogram(formula= 	dataStnMonDiffSpdf@data[,1] ~ 1, 
		input_data= dataStnMonDiffSpdf,
		verbose=	TRUE)
plot(dataStnMonDiffGstatVgm, 
		model=dataStnMonDiffAutoVgmModFit$var_model, 
		main="[automap autofitVariogram] Variogram: (Station-Satellite) Differences",
		sub=paste("Model:", 		dataStnMonDiffAutoVgmModFit$var_model[2,1], 
				",   Nugget: ", 	round(dataStnMonDiffAutoVgmModFit$var_model[1,2]), 
				",   Full Sill: ", 	round(dataStnMonDiffAutoVgmModFit$var_model[2,2] + dataStnMonDiffAutoVgmModFit$var_model[1,2]), 
				",   Range: ", 		round(dataStnMonDiffAutoVgmModFit$var_model[2,3], digits=2), 
				sep=""),
		col="black")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_automap_fit_Rain_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)


# Use automap parameters for Exponential model to start fit.variogram from gstat

dataStnMonStnAutoGstatVgmModFit	= fit.variogram(dataStnMonStnGstatVgm,  
												vgm(dataStnMonStnAutoVgmModFit$var_model[2,2],
													"Exp",
													dataStnMonStnAutoVgmModFit$var_model[2,3],
													dataStnMonStnAutoVgmModFit$var_model[1,2]))
plot(dataStnMonStnGstatVgm, model=dataStnMonStnAutoGstatVgmModFit, col="black", main="[automap & gstat] Variogram: Stations")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_auto_gstat_fit_Rain_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

dataStnMonSatAutoGstatVgmModFit	= fit.variogram(dataStnMonSatGstatVgm,  
												vgm(dataStnMonSatAutoVgmModFit$var_model[2,2],
														"Exp",
														dataStnMonSatAutoVgmModFit$var_model[2,3],
														dataStnMonSatAutoVgmModFit$var_model[1,2]))
plot(dataStnMonSatGstatVgm, model=dataStnMonSatAutoGstatVgmModFit, col="black", main="[automap & gstat] Variogram: Satellite (FTIP)")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_auto_gstat_fit_Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

dataStnMonDiffAutoGstatVgmModFit	= fit.variogram(dataStnMonDiffGstatVgm,  
											vgm(dataStnMonDiffAutoVgmModFit$var_model[2,2],
													"Exp",
													dataStnMonDiffAutoVgmModFit$var_model[2,3],
													dataStnMonDiffAutoVgmModFit$var_model[1,2]))
plot(dataStnMonDiffGstatVgm, model=dataStnMonDiffAutoGstatVgmModFit, col="black", main="[automap & gstat] Variogram: (Station - Satellite) Differences")
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_auto_gstat_fit_Rain_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)








# Use automap parameters for Exponential model to start variog from geoR

dataStnMonStnAutoGeoVgmModFit	= variofit(dataStnMonGeoVgmStn, 
											ini.cov.pars=c(dataStnMonStnAutoVgmModFit$var_model[2,2],
															dataStnMonStnAutoVgmModFit$var_model[2,3]))
plot(dataStnMonGeoVgmStn, main="[automap & geoR] Variogram: Stations")
lines(dataStnMonStnAutoGeoVgmModFit)
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_auto_geoR_fit_Rain_Stn_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
											

dataStnMonSatAutoGeoVgmModFit	= variofit(dataStnMonGeoVgmSat, 
		ini.cov.pars=c(dataStnMonSatAutoVgmModFit$var_model[2,2],
				dataStnMonSatAutoVgmModFit$var_model[2,3]))
plot(dataStnMonGeoVgmSat, main="[automap & geoR] Variogram: Satellite (FTIP)")
lines(dataStnMonSatAutoGeoVgmModFit)
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_auto_geoR_fit_Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)


dataStnMonDiffAutoGeoVgmModFit	= variofit(dataStnMonGeoVgmDiff, 
		ini.cov.pars=c(dataStnMonDiffAutoVgmModFit$var_model[2,2],
				dataStnMonDiffAutoVgmModFit$var_model[2,3]))
plot(dataStnMonGeoVgmDiff, main="[automap & geoR] Variogram: (Station - Satellite) Differences")
lines(dataStnMonDiffAutoGeoVgmModFit)
dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Vgm_auto_geoR_fit_Rain_Diff_", sprintf("%02d", month),".png", sep=""), width=480, height=480)


# Now
# load the FTIP as a raster, 
# compute a variogram in a moving window (2 degrees), 
# either in adjacent or overlapping windows (1 degree),
# collect the nugget, range, and sill parameters, 
# map the parameters 



#	movingFun(x, n, fun=mean, type='around', circular=FALSE, na.rm=FALSE) 
#	focal(x, w=3, fun, filename='', na.rm=FALSE, pad=FALSE, padValue=NA, NAonly=FALSE, ...) 
#	getValuesBlock(x, row, ...)
#	getValuesFocal(x, row, nrows, ngb, ...)
#	gridDistance(x, origin, omit=NULL, filename="", ...) 
#	pairs(x, ...)
#	plot(x, col=rev(terrain.colors(25)), maxpixels=100000, axes = TRUE, xlab="", ylab="", alpha=1, useRaster=TRUE, ...) 
#	pointDistance(p1, p2, longlat, ...)
#	
#	canProcessInMemory(x, n=4)
#	closeConnection(x)
#	openConnection(x, silent=FALSE)
#	pbCreate(nsteps, progress, style=3, ...)
#	pbStep(pb, step=NULL, label='')
#	pbClose(pb, timer)
#	getCluster()
#	returnCluster()
#	
#	rasterToPoints(x, fun=NULL, spatial=FALSE, ...)
#	
#	replace:
#	x[i] <- value 
#	x[i,j] <- value 
#	
#	scalebar(d, xy = NULL, type = "line", divs = 2, below = "", lonlat = NULL, label, adj=c(0.5, -0.5), lwd = 2, ...)
#	setValues(x, values, layer) 
#	values(x) <- value
#	spplot(obj, ..., maxpixels=50000, as.table=TRUE) 
#	
#	rasterVis package 
#	
#	subs(x, y, by=1, which=2, subsWithNA=TRUE, filename='', ...)
#	writeRaster(x, filename, format, ...)
#	
#	writeStart(x, filename, ...)
#	writeValues(x, v, start)
#	writeStop(x)
#	
#	xFromCol(object, colnr)
#	yFromRow(object, rownr)
#	xyFromCell(object, cell, spatial=FALSE)
#	xFromCell(object, cell)
#	yFromCell(object, cell)
#	
#	zonal(x, zones, stat='mean', digits=0, na.rm=TRUE, progress='') 
#	
#	setZ(x, z, name='time')
#	getZ(x)
#
#	vgram.matrix {fields}

#	If the data object passed to gstat:variogram is a spatial object (your data is aSpatialGridDataFrame) then you do not need to specify the locations, as these are contained in the data.
#	
#	However, clearly the NA values are causing a problem, so if we force the grid object to be a SpatialPointsDataFrame, this will remove the NA values
#	
#	im contains the data https://gist.github.com/2780792
#	
#	library(gstat)
#	point_data <- as(im, 'SpatialPointsDataFrame')
#	gstat_variogram <- variogram(band1 ~ 1, data = point_data)
#	
#	To use geoR
#	
#	library(geoR)
#	geor_variogram <- variog(coords = coordinates(point_data), 
#			data = point_data@data$band1)
#	
#	or even more simply (as geoR works with objects of class geodata and contains the function as.geodata to convert from a SpatialPointsDataFrame to geodata object
#	
#	geor_variogram <- variog(as.geodata(point_data))

#	r <- raster(ncols=100, nrows=100)
#	cells <- c(3:5, 210)
#	r <- rasterFromCells(r, cells)
#	cbind(1:ncell(r), getValues(r))





#	inputRasterImg			=	FTIPMonImg

rasterMovWinVgmPar	=	function (	inputRasterImg,
									outputRasterNuggetName, outputRasterRangeName, outputRasterSillName,
									cellLatDim, cellLonDim, 
									winLatDimPre, winLatDimPost,
									winLonDimPre, winLonDimPost)
	{

	# Create output rasters for
	# Nugget, Range, and Sill values
	# initialized as NA values.
	outputRasterNugget	= 	inputRasterImg
	outputRasterNugget	=	setValues(outputRasterNugget, values=matrix(NA, ncol=ncell(outputRasterNugget)))
	outputRasterSill	= 	outputRasterNugget
	outputRasterRange	= 	outputRasterNugget
	
#	# Coerce inputRaster to 'SpatialPointsDataFrame'
#	inputSpdf			=	rasterToPoints(inputRaster, fun=NULL, spatial=TRUE)
#	
#	# Coerce inputSpdf to 'geodata'
#	inputGeo 			=	as.geodata(inputSpdf, data.col = 1)
	
	
#	winLatDimPre	=	WinLatDimPre
#	winLatDimPost	=	WinLatDimPost
#	winLonDimPre	=	WinLonDimPre
#	winLonDimPost	=	WinLonDimPost
	
#	cellLatDim	=	CellLatDim
#	cellLonDim	=	CellLonDim
	
#	winRowDimPre	=	WinLatDimPre  / CellLatDim
#	winRowDimPost	=	WinLatDimPost / CellLatDim
#	winColDimPre	=	WinLonDimPre  / CellLonDim
#	winColDimPost	=	WinLonDimPost / CellLonDim	
	
	winRowDimPre	=	winLatDimPre  / cellLatDim
	winRowDimPost	=	winLatDimPost / cellLatDim
	winColDimPre	=	winLonDimPre  / cellLonDim
	winColDimPost	=	winLonDimPost / cellLonDim	
	
	# For loops
	for (row in (1+winRowDimPre):(nrow(inputRasterImg)-winRowDimPost))			# row=1+winRowDimPre : row=nrow(inputRasterImg)-winRowDimPost
		{
		for (col in (1+winColDimPre):(ncol(inputRasterImg)-winColDimPost))		# col=1+winColDimPre : col=ncol(inputRasterImg)-winColDimPost
			{
			
			cat("============================\n")
			cat("Starting at new location ...\n")
			cat(paste("Row :", row, "in rows", 1+winRowDimPre, "to", nrow(inputRasterImg)-winRowDimPost, "\n"))
			cat(paste("Col :", col, "in cols", 1+winColDimPre, "to", ncol(inputRasterImg)-winColDimPost, "\n"))
		
#			# Subset window from each point data set (SPDF and geodata)
#			inputSpdfWin	=	inputSpdf	[row-winRowDimPre:row+winRowDimPost, col-winColDimPre:col+winColDimPost]
#			inputGeoWin		=	inputGeo	[row-winRowDimPre:row+winRowDimPost, col-winColDimPre:col+winColDimPost]
#print (paste(row-winRowDimPre, row+winRowDimPost, col-winColDimPre, col+winColDimPost))
			winExtent		=	extent	(inputRasterImg, row-winRowDimPre, row+winRowDimPost, col-winColDimPre, col+winColDimPost)
			inputRasterWin	=	crop	(inputRasterImg, winExtent)

			# Coerce inputRasterWin to 'SpatialPointsDataFrame'
			inputSpdfWin	=	rasterToPoints(inputRasterWin, fun=NULL, spatial=TRUE)
			
			# Coerce inputSpdfWin to 'geodata'
			inputGeoWin 	=	as.geodata(inputSpdfWin, data.col = 1)			
			
			# Calculate variogram parameters for values of subset window (using 'automap' package, 'autofitVariogram' function)
			autoVgmPar	= autofitVariogram(	formula		=	inputSpdfWin@data[,1] ~ 1, 
											input_data	=	inputSpdfWin,
											verbose		=	TRUE)
			
			
			# Calculate variogram for values of subset window (using 'geoR' package, 'variog' function)
			geoVgm		= variog(	inputGeoWin, 
									coords	=	inputGeoWin$coords, 
									data	=	inputGeoWin$data, 
									max.dist=	(1/2) * max((winLatDimPre+winLatDimPost),(winLonDimPre+winLonDimPost)))
			
									
			# Calculate variogram parameters for values of subset window (using 'geoR' package, 'variofit' function)
			fitVgmPar	= variofit(	geoVgm, 
									ini.cov.pars	=	c(	autoVgmPar$var_model[2,2],
															autoVgmPar$var_model[2,3]))
												
			valueNugget	= fitVgmPar$nugget
			valueSill	= fitVgmPar$cov.pars[1] + fitVgmPar$nugget
			valueRange	= fitVgmPar$cov.pars[2]
						
			outputRasterNugget	[row,col]	=	valueNugget
			outputRasterSill	[row,col]	=	valueSill
			outputRasterRange	[row,col]	=	valueRange

			assign(deparse(substitute(outputRasterNuggetName)), outputRasterNugget, envir = .GlobalEnv)
			assign(deparse(substitute(outputRasterSillName)), 	outputRasterSill, envir = .GlobalEnv)
			assign(deparse(substitute(outputRasterRangeName)), 	outputRasterRange, envir = .GlobalEnv)
			
			plot	(FTIPMonImg,
					col=colorRampPalette(c("black", "white"))(25),
					main= "Rainfall: Satellite (FTIP)",
					xlab= "Longitude", ylab= "Latitude")
			map		('worldHires', add=TRUE, col= "green")
			dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Rain_Sat_", sprintf("%02d", month),"_working.png", sep=""), width=480, height=480)
			
			plot	(outputRasterNugget,
					col=colorRampPalette(c("black", "white"))(25),
					main= "Rainfall:  Mov. Win. Variogram Nugget",
					xlab= "Longitude", ylab= "Latitude")
			map		('worldHires', add=TRUE, col= "green")
			dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Nugget_", sprintf("%02d", month),"_working.png", sep=""), width=480, height=480)
			
			plot	(outputRasterSill,
					col=colorRampPalette(c("black", "white"))(25),
					main= "Rainfall: Mov. Win. Variogram Sill",
					xlab= "Longitude", ylab= "Latitude")
			map		('worldHires', add=TRUE, col= "green")
			dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Sill_", sprintf("%02d", month),"_working.png", sep=""), width=480, height=480)
			
			plot	(outputRasterRange,
					col=colorRampPalette(c("black", "white"))(25),
					main= "Rainfall: Mov. Win. Variogram Range",
					xlab= "Longitude", ylab= "Latitude")
			map		('worldHires', add=TRUE, col= "green")
			dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Range_", sprintf("%02d", month),"_working.png", sep=""), width=480, height=480)
			
			}
		}
		
#	assign(deparse(substitute(outputRasterNuggetName)), outputRasterNugget, envir = .GlobalEnv)
#	assign(deparse(substitute(outputRasterSillName)), 	outputRasterSill, envir = .GlobalEnv)
#	assign(deparse(substitute(outputRasterRangeName)), 	outputRasterRange, envir = .GlobalEnv)

	# Write output raster results
	writeRaster	(FTIPMonSmallTest, 		paste(DirPrecipVgmPar, outputDir, "\\", "Raster_Rain_Small_Test_", 		sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
	writeRaster	(FTIPMonImg, 			paste(DirPrecipVgmPar, outputDir, "\\", "Raster_Rain_Sat_", 			sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
	writeRaster	(FTIPMonNugget, 		paste(DirPrecipVgmPar, outputDir, "\\", "Raster_Vgm_Nugget_", 			sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
	writeRaster	(FTIPMonSill, 			paste(DirPrecipVgmPar, outputDir, "\\", "Raster_Vgm_Sill_", 			sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
	writeRaster	(FTIPMonRange, 			paste(DirPrecipVgmPar, outputDir, "\\", "Raster_Vgm_Range_", 			sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
	writeRaster	(FTIPMonNugget^(1/5),	paste(DirPrecipVgmPar, outputDir, "\\", "Raster_Vgm_Nugget_Trnsfm_",	sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
	writeRaster	(FTIPMonSill^(1/5), 	paste(DirPrecipVgmPar, outputDir, "\\", "Raster_Vgm_Sill_Trnsfm_", 		sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
	writeRaster	(FTIPMonRange^(1/5),	paste(DirPrecipVgmPar, outputDir, "\\", "Raster_Vgm_Range_Trnsfm_", 	sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
	
	
	
	# Plot the input data
	plot	(FTIPMonSmallTest,
			col=colorRampPalette(c("black", "white"))(25),
			main= "Rainfall: Satellite (FTIP)",
			xlab= "Longitude", ylab= "Latitude")
	map		('worldHires', add=TRUE, col= "green")
	dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Rain_Small_Test_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
	
	plot	(FTIPMonImg,
			col=colorRampPalette(c("black", "white"))(25),
			main= "Rainfall: Satellite (FTIP)",
			xlab= "Longitude", ylab= "Latitude")
	map		('worldHires', add=TRUE, col= "green")
	dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Rain_Sat_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
	
	# Plot variogram parameters
	plot	(FTIPMonNugget,
			col=colorRampPalette(c("black", "white"))(25),
			main= "Rainfall: Mov. Win. Variogram Nugget",
			xlab= "Longitude", ylab= "Latitude")
	map		('worldHires', add=TRUE, col= "green")
	dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Nugget_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
	
	plot	(FTIPMonSill,
			col=colorRampPalette(c("black", "white"))(25),
			main= "Rainfall: Mov. Win. Variogram Sill",
			xlab= "Longitude", ylab= "Latitude")
	map		('worldHires', add=TRUE, col= "green")
	dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Sill_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
	
	plot	(FTIPMonRange,
			col=colorRampPalette(c("black", "white"))(25),
			main= "Rainfall: Mov. Win. Variogram Range",
			xlab= "Longitude", ylab= "Latitude")
	map		('worldHires', add=TRUE, col= "green")
	dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Range_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
	
	
	# Plot variogram parameters, transformed (^1/10)
	plot	(FTIPMonNugget^(1/10),
			col=colorRampPalette(c("black", "white"))(25),
			main= "Rainfall: Mov. Win. Variogram Nugget (^1/10)",
			xlab= "Longitude", ylab= "Latitude")
	map		('worldHires', add=TRUE, col= "green")
	dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Nugget_Trnsfm_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
	
	plot	(FTIPMonSill^(1/10),
			col=colorRampPalette(c("black", "white"))(25),
			main= "Rainfall: Mov. Win. Variogram Sill (^1/10)",
			xlab= "Longitude", ylab= "Latitude")
	map		('worldHires', add=TRUE, col= "green")
	dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Sill_Trnsfm_", sprintf("%02d", month),".png", sep=""), width=480, height=480)
	
	plot	(FTIPMonRange^(1/10),
			col=colorRampPalette(c("black", "white"))(25),
			main= "Rainfall: Mov. Win. Variogram Range (^1/10)",
			xlab= "Longitude", ylab= "Latitude")
	map		('worldHires', add=TRUE, col= "green")
	dev.print(png,file=paste(DirPrecipVgmPar, outputDir, "Raster_Vgm_Range_Trnsfm_", sprintf("%02d", month),".png", sep=""), width=480, height=480)

	}
	
	

#	plot	(FTIPMonImg,
##		ext=FTIPMonImg@extent,
##		col=rev(rainbow(25)),
#			col=colorRampPalette(c("black", "white"))(25),
##		col=colorRampPalette(c("blue", "white", "red"))(25),
##		col=colorRampPalette(c("brown","gray80","green"))(25)[cut(sqrt(FTIPMonImg@data@values),25)],
#			main= "Rainfall: Satellite (FTIP)",
#			xlab= "Longitude", ylab= "Latitude")
#	map		('worldHires', add=TRUE, col= "green")
	
	
#levelplot	(FTIPMonImg,
#			margin=FALSE,
#			col.regions=terrain.colors(25),
#			cuts=25,
#			main= "Rainfall: Satellite (FTIP)",
#			xlab= "Longitude", ylab= "Latitude")
#map			('worldHires', add=TRUE, col= "green")
	
	
#rm(FTIPMonNugget, FTIPMonRange, FTIPMonSill)

# Run moving window function
# on testing data (5 degree X 5 degree; minus 1 degree on each edge -> 3 dg. X 3 dg.)
# Took about 4.5 hours for 9 square degrees, or 0.5 hours for each sq. degree
# For 360x100 = 36000 sq degrees, that would be 18000 hours in serial
# or 18000 hours / 80 hours = 225 parallel processes.
# Need to reduce automap parameter estimation, remove geoR fitting, and/or restrict to land areas
rasterMovWinVgmPar	(	FTIPMonSmallTest,
						FTIPMonNugget, FTIPMonRange, FTIPMonSill,
						CellLatDim, CellLonDim, 
						WinLatDimPre, WinLatDimPost,
						WinLonDimPre, WinLonDimPost)

# Run moving window function
# on CACOL data (24 degree X 25 degree)
rasterMovWinVgmPar	(	FTIPMonImg,
						FTIPMonNugget, FTIPMonRange, FTIPMonSill,
						CellLatDim, CellLonDim, 
						WinLatDimPre, WinLatDimPost,
						WinLonDimPre, WinLonDimPost)

















##	End of for loops for new CACOL variogram code.		
#		}
#		
#}
	















# Global latitude calculations for 'defineRegions' function
LatRangeNominal	= LatMaxNominal - LatMinNominal
LatRegionCount	= LatRangeNominal/LatRegionDim
LatRegionBreaksNominal = seq(LatMinNominal, LatMaxNominal, LatRegionDim)

# Global longitude calculations for 'defineRegions' function
LonRangeNominal	= LonMaxNominal - LonMinNominal
LonRegionCount	= LonRangeNominal/LonRegionDim
LonRegionBreaksNominal = seq(LonMinNominal, LonMaxNominal, LonRegionDim)

# The following code is wrapped in a function to use try() to handle errors in 'Run_vario.R'
# in case the variogram parameter file hasn't yet been generated with 'Run_vario.R'.
# These preliminary data calculations are for 'Run_krige.R',
# particularly for the mosaic and mask operations.
# The variogram parameter file should be created before kriging.

vgmParamList <- function () {
	# Read variogram parameters for systematic collection of
	# FileTagRegions and DirTagMonthsKrg (see below) for use
	# in mosaic (krgMosaic, run_krgMosaic)
	# and mask (krgLandMask, run_krgLandMask) functions
	DataVgmPar = 	read.csv(paste(DirPrecipVgmPar, FilePrecipVgmPar, sep=""),
					header=TRUE, stringsAsFactors=FALSE, strip.white=TRUE)
	
	
	# For specifying regions
	# in mosaic (krgMosaic, run_krgMosaic)
	# and mask (krgLandMask, run_krgLandMask) functions
	#
	## Hard-coded specification; confirm that there is no comma after last entry
	#		FileTagRegions = c(
	#		"Region_0_50N_180W_120W",
	#		"Region_0_50N_120W_60W",
	#		"Region_0_50N_60W_0",
	#		"Region_0_50N_0_60E",
	#		"Region_0_50N_60E_120E",
	#		"Region_0_50N_120E_180E",
	#		"Region_50S_0_180W_120W",
	#		"Region_50S_0_120W_60W",
	#		"Region_50S_0_60W_0",
	#		"Region_50S_0_0_60E",
	#		"Region_50S_0_60E_120E",
	#		"Region_50S_0_120E_180E"
	#		)
	#
	## Systematic specification from variogram parameters file
	fileTagRegions = unique(DataVgmPar$Region)
	assign("FileTagRegions", fileTagRegions, envir = .GlobalEnv)
	
	
	# For specifying months
	# in mosaic (krgMosaic, run_krgMosaic)
	# and mask (krgLandMask, run_krgLandMask) functions
	#
	## Hard-coded specification; confirm no comma after last entry
	#		DirTagMonthsKrg = c(1,2,3,4,5,6,7,8,9,10,11,12)
	#
	## Systematic specification from variogram parameters file
	dirTagMonthsKrg = unique(DataVgmPar$Month)
	assign("DirTagMonthsKrg", dirTagMonthsKrg, envir = .GlobalEnv)
	
	# For specifying kriged years
	# in mosaic (krgMosaic, run_krgMosaic)
	# and mask (krgLandMask, run_krgLandMask) functions
	fileTagYearsKrg	= seq.int(YearBegKrg, YearEndKrg)
	assign("FileTagYearsKrg", fileTagYearsKrg, envir = .GlobalEnv)
	
	# For specifying the kriging output types
	# in mosaic (krgMosaic, run_krgMosaic)
	# and mask (krgLandMask, run_krgLandMask) functions
	fileTagKrgOutputPrecip	= c(FilePrefixPrecipKrgEstMm, FilePrefixPrecipKrgEstPc, FilePrefixPrecipKrgErrMm, FilePrefixPrecipKrgErrPc)
	assign("FileTagKrgOutputPrecip", fileTagKrgOutputPrecip, envir = .GlobalEnv)
	
	fileTagKrgOutputTemp	= c(FilePrefixTempKrgEstDg, FilePrefixTempKrgErrDg)
	assign("FileTagKrgOutputTemp", fileTagKrgOutputTemp, envir = .GlobalEnv)
}

#try(vgmParamList(), silent=TRUE)






########################################################################
########################################################################
########################################################################
# CODE FOR FITTING EMPIRICAL VARIOGRAM MODELS
# Some VARIOGRAM PREPARATION is in KRIGING section (e.g. prepStnAnmly).
########################################################################

# Define function "fit_vario"
# region = region of interest
# mmm = month definition (i.e. jan, feb, mar, apr, may, jun, jul, aug, sep, oct, nov, dec)
# mode = "t" for temperature or "p" for precipitation
fit_vario = function(mode,region,month) {
	
	# Source file 'vario_krige.R' has required variables and functions (e.g. avgVgm)
	source(paste(DirPrecipSourceVgmKrg, FilePrecipSourceVgmKrg, sep=""))
	
	# Maybe source defineRegions from external file
	# (e.g. defineRegions.R), with no duplicate in vario_krige.R,
	# or run ONCE in Run_vario.R (The latter implementation is currently implemented.)
#	defineRegions()
	
	cat(paste("Mode: ", mode, "\n", sep=""))
	cat(paste("Region: ", region$regionName, "\n", sep="")) 
	cat(paste("Month: ", month, "\n", sep=""))
	
	# Prepare station data
	prepData = prepStnAnmly (mode,region,month)
	
	locs = cbind(prepData$meta$Lon,prepData$meta$Lat)
	colnames(locs) = c("lon","lat")
	
	# Default starting variogram model parameters	
	if( mode == "p" ) { 			# p is for precipitation
		range  = 6
		nugget = .01
		sill   = 0.032
	}
	
	if( mode == "t" ) { 			# t is for temperature
		range  = 6
		nugget = 0.1
		sill   = 0.5
	}
	
	
	while(1) {
		
		# create variogram
		if( mode == "p" ) vgr=avgVgm(prepData$pc^0.5,locs,region,YearBegVgm,YearEndVgm,range,nugget,sill)
		if( mode == "t" ) vgr=avgVgm(prepData$pc    ,locs,region,YearBegVgm,YearEndVgm,range,nugget,sill)		
		
		# Option A: To just return the variogram model, use the next line:
		return(vgr)
		
		# Option B: To re-run the whole (time-consuming) process,
		#			use this last section of code:
#		print("Rerun with new paramaters? (y/n)")
#		ok = readline()
#		if( ok == "y" ) {
#			print("Enter Range:")
#			range = as.numeric(readline())
#			print("Enter Nugget:")
#			nugget = as.numeric(readline())
#			print("Enter Sill")
#			sill = as.numeric(readline())
#			print(c("Range:",range))
#			print(c("Nugget:",nugget))
#			print(c("Sill:",sill))
#		}
#		if( ok == "n" ) {
#			return(vgr)
#		}
		
		
	}
}



#
# avgVgm
# 'data' is matrix of (station) values with each column defining a different (year) time.
#
avgVgm = function(data,locs,region,begYear,endYear,range,nugget,sill) {
	
	require(geoR)	
	require(gstat)
	
	data = data.frame(data)						# coerce to data frame for ease of using column names
	colnames(data) = DataFileColNames	
	colYearBegVgm = which(names(data)==begYear)	# find column index for beginning year of variogram
	colYearEndVgm = which(names(data)==endYear)	# find column index for end year of variogram	
	
	lon = as.numeric(locs[,1])
	lat = as.numeric(locs[,2])
	obs = as.numeric(data[,colYearBegVgm])
	
	n   = length(lon)   
	tmp = matrix(NA,nrow=n,ncol=3)
	
	tmp[,1] = lon
	tmp[,2] = lat
	tmp[,3] = obs
	
	gid = !is.na(tmp[,3]) 
	
	gdat    = as.geodata(dedup(tmp[gid,]))
	colnames(tmp) = c("lon","lat","obs")
	
	gtmp = data.frame(dedup(tmp[gid,]))
	# Alternative variogram model below for other purposes.
	# bing = variogram(obs~lon+lat,loc=~lon+lat,gtmp,cutoff=8)
	bing = variogram(obs~1,loc=~lon+lat,gtmp,cutoff=10)
	vg = bing$gamma								
	ug = bing$dist
	vg[] = 0								
	ug[] = 0		
#	plot(bing$dist,bing$gamma)	# Supress plot for speed
	
	cnt = 0	
	for( i in colYearBegVgm:colYearEndVgm ) {
		
		tmp[,3] = as.numeric(data[,i])
		
		gid = !is.na(tmp[,3]) & tmp[,3] <= 2000
		gdat    = as.geodata(dedup(tmp[gid,]))
		gtmp  = data.frame(dedup(tmp[gid,]))
		# Alternative variogram model below for other purposes.
		# bing2 = variogram(obs~lon+lat,loc=~lon+lat,gtmp,cutoff=6)
		bing2 = variogram(obs~1,loc=~lon+lat,gtmp,cutoff=10)
		vg = vg + bing2$gamma					
		ug = ug + bing2$dist					
#		plot(bing2$dist,bing2$gamma,main=paste("Variogram for",yr))	# Supress plot for speed
		cnt = cnt + 1							
		cat(colnames(data)[i])
		cat("\n")
	}
	
	
	bing$dist  = ug / cnt
	bing$gamma = vg / cnt
	
	plot(bing$dist,bing$gamma,main="Variogram based on average semivariances",ylim=c(0,max(bing$gamma))	)
	
	
	vmod = vgm(sill-nugget,VgmModelForm,range,nugget)
	
	lines(variogramLine(vmod,MaxVarioDist))
	return(vmod)
}





########################################################################
########################################################################
########################################################################
# CODE FOR KRIGING, MOSAICS, AND MASKING.
########################################################################
# Define function "run_run_mons_p" to
# 1) set the region of interest (region)
# 2) run the "run_mon" function several times
#    using the region of interest (region), the months (mm = 01 through 12), 
#    and the *** empirically modeled variogram parameters (from Run_vario.R) ***
#    range, nugget, sill values
#    >run_mon = function(region,mm,range,nugget,sill)

run_run_mon_p = function() {
	
	# Source file 'vario_krige.r' has other required functions
	source(paste(DirPrecipSourceVgmKrg, FilePrecipSourceVgmKrg, sep=""))
	
	# Define regions based upon global variables
	defineRegions()
	
	# Read variogram parameters for use in kriging
	dataVgmPar = 	read.csv(paste(DirPrecipVgmPar, FilePrecipVgmPar, sep=""),
			header=TRUE, stringsAsFactors=FALSE, strip.white=TRUE)
	
	# Run kriging per month iteratively per variogram parameters
	for(i in 1:nrow(dataVgmPar)) {
		
		run_mon(dataVgmPar$Mode[i],
				get(dataVgmPar$Region[i]),
				dataVgmPar$Month[i],
				dataVgmPar$Nugget[i],
				dataVgmPar$Sill[i],
				dataVgmPar$Range[i]
		)
	}
}



#
# This function creates a control structure defining our interpolation region
#
defineRegions = function () {
	for (regionLat in 1:LatRegionCount) {
		for (regionLon in 1:LonRegionCount) {
			
			regionLatMinNominal = LatRegionBreaksNominal[regionLat]
			regionLatMaxNominal = LatRegionBreaksNominal[regionLat+1]
			regionLonMinNominal = LonRegionBreaksNominal[regionLon]
			regionLonMaxNominal = LonRegionBreaksNominal[regionLon+1]
			
			if (regionLatMinNominal <  0) {regionLatMinHemi = "S"}
			if (regionLatMinNominal == 0) {regionLatMinHemi = ""}
			if (regionLatMinNominal >  0) {regionLatMinHemi = "N"}
			
			if (regionLatMaxNominal <  0) {regionLatMaxHemi = "S"}
			if (regionLatMaxNominal == 0) {regionLatMaxHemi = ""}
			if (regionLatMaxNominal >  0) {regionLatMaxHemi = "N"}
			
			if (regionLonMinNominal <  0) {regionLonMinHemi = "W"}
			if (regionLonMinNominal == 0) {regionLonMinHemi = ""}
			if (regionLonMinNominal >  0) {regionLonMinHemi = "E"}
			
			if (regionLonMaxNominal <  0) {regionLonMaxHemi = "W"}
			if (regionLonMaxNominal == 0) {regionLonMaxHemi = ""}
			if (regionLonMaxNominal >  0) {regionLonMaxHemi = "E"}
			
			regionName=paste("Region_",
					abs(regionLatMinNominal), regionLatMinHemi,"_",
					abs(regionLatMaxNominal), regionLatMaxHemi,"_",
					abs(regionLonMinNominal), regionLonMinHemi,"_",
					abs(regionLonMaxNominal), regionLonMaxHemi,
					sep="")
			
			
			regionLatMin			= regionLatMinNominal + (CellLatDim/2)
			regionLatMax			= regionLatMaxNominal - (CellLatDim/2)
			regionLatCellCntrs		= seq(regionLatMax, regionLatMin, -CellLatDim) # max to min, to prepare for listing from upper-left to lower-right
			regionLatCellCntrsCnt	= length(regionLatCellCntrs)
			
			regionLonMin			= regionLonMinNominal + (CellLonDim/2)
			regionLonMax			= regionLonMaxNominal - (CellLonDim/2)
			regionLonCellCntrs		= seq(regionLonMin, regionLonMax, CellLonDim) # min to max, to prepare for listing from upper-left to lower-right
			regionLonCellCntrsCnt	= length(regionLonCellCntrs)
			
			regionData = list(
					
					regionName				= regionName,
					
					regionLatMinNominal		= regionLatMinNominal,
					regionLatMaxNominal		= regionLatMaxNominal,
					regionLatMin			= regionLatMin,
					regionLatMax			= regionLatMax,
					regionLatCellDim		= CellLatDim,
					regionLatCellCntrs		= regionLatCellCntrs,
					regionLatCellCntrsCnt	= regionLatCellCntrsCnt,
					
					regionLonMinNominal		= regionLonMinNominal,
					regionLonMaxNominal		= regionLonMaxNominal,				
					regionLonMin			= regionLonMin,
					regionLonMax			= regionLonMax,			
					regionLonCellDim		= CellLonDim,								
					regionLonCellCntrs		= regionLonCellCntrs,
					regionLonCellCntrsCnt	= regionLonCellCntrsCnt
			
			)	
			
			assign(regionName, regionData, pos = .GlobalEnv)			
		}
	}
}



########################################################################
# Define function "run_mon" to run kriging per months
# mode		= variable mode ("p" for precip, "t" for temp)
# region    = region of interest
# month     = month (in two digits as text, i.e. "01" through "12")
# nugget  	= nugget of empirically modeled variogram 
# sill    	= sill of empirically modeled variogram 
# range   	= range of empirically modeled variogram 

run_mon = function(mode,region,month,nugget,sill,range) {
	
	# Load the package "gstat"
	require(gstat)
	
	# Source file 'vario_krige.r' has other required functions
	source(paste(DirPrecipSourceVgmKrg, FilePrecipSourceVgmKrg, sep=""))
	
	# Print function parameters
	cat(paste("Mode: ", 				mode, 				"\n", sep=""))
	cat(paste("Region: ", 				region$regionName, 	"\n", sep="")) 
	cat(paste("Month: ", 				month, 				"\n", sep=""))
	cat(paste("Variogram nugget: ", 	nugget, 			"\n", sep=""))
	cat(paste("Variogram sill: ", 		sill, 				"\n", sep=""))
	cat(paste("Variogram range: ", 		range, 				"\n", sep=""))
	
	# Create variogram model.
	# Note: partial sill needed, instead of sill, hence, (full)sill-nugget.
	# vgm(psill, model, range, nugget, add.to, anis, kappa = 0.5, ..., covtable)
	vgmModel = vgm(sill-nugget,VgmModelForm,range,nugget)
	
	# Prepare data for variogram fitting and kriging operations.
	prepData = prepStnAnmly (mode,region,month)
	
	# Prepare background reference data to guide interpolation ...
	
	# Create a set of lat/lon values to interpolate to
	toLonLat = regionLocs(region)
	
	# Create output directory from global variable and month
	if( mode == "p" ) {
		outputDir = sprintf("%02d", month)
		dir.create(file.path(DirPrecipKrige,outputDir), recursive=TRUE, showWarnings=FALSE)
	}
	
	if( mode == "t" ) {
		outputDir = sprintf("%02d", month)
		dir.create(file.path(DirTempKrige,outputDir), recursive=TRUE, showWarnings=FALSE)
	}
	
	# Create a set of lat/lon values to interpolate from
	fromLonLat = cbind(prepData$meta$Lon,prepData$meta$Lat)
	colnames(fromLonLat) = c("lon","lat")
	
	# Read background reference data
	RefData = readRefData(month) 
	
	# Identify years of input data to krige	
	colYearBegKrg = which(colnames(prepData$pc)==YearBegKrg)
	colYearEndKrg = which(colnames(prepData$pc)==YearEndKrg)
	
	# Krige iteratively per year
	for( i in colYearBegKrg:colYearEndKrg ) {
		
		# The current year that's being kriged
		yearKrg = (YearBegData-1)+i
		cat(paste("Kriging year: ", yearKrg, "\n", sep=""))
		
		if( mode == "p" ) {
			
			# Block kriging of current month, current year, per variogram paramters and prepared data.
			kr = krigblock(region,vgmModel,toLonLat,fromLonLat,prepData$pc[,i]^0.5,c(CellLonDim,CellLatDim)) 
			pr.pc = kr$predict
			pr.er = kr$krige.var^0.5
			
			# Convert kriging from (square root) percent to millimeters
			krmm  = pct2mm(region,RefData,pr.pc,pr.er,toLonLat)
			
			# Output kriged prediction values in millimeters
			krmm_mmkr = matrix (krmm$mmkr, region$regionLatCellCntrsCnt, region$regionLonCellCntrsCnt, byrow=TRUE)
			krmm_mmkr = raster (krmm_mmkr, xmn=region$regionLonMinNominal, xmx=region$regionLonMaxNominal, ymn=region$regionLatMinNominal, ymx=region$regionLatMaxNominal, crs="+proj=lonlat +datum=WGS84")
			writeRaster	(krmm_mmkr, paste(DirPrecipKrige, outputDir, "\\", FilePrefixPrecipKrgEstMm, "_", region$regionName, "_", yearKrg, "_", sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
			
			# Output kriged standard error values in millimeters
			krmm_mmse = matrix (krmm$mmse, region$regionLatCellCntrsCnt, region$regionLonCellCntrsCnt, byrow=TRUE)
			krmm_mmse = raster (krmm_mmse, xmn=region$regionLonMinNominal, xmx=region$regionLonMaxNominal, ymn=region$regionLatMinNominal, ymx=region$regionLatMaxNominal, crs="+proj=lonlat +datum=WGS84")
			writeRaster	(krmm_mmse, paste(DirPrecipKrige, outputDir, "\\", FilePrefixPrecipKrgErrMm, "_", region$regionName, "_", yearKrg, "_", sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
			
			# Output kriged prediction values in percent
			krmm_pckr = matrix (krmm$pckr, region$regionLatCellCntrsCnt, region$regionLonCellCntrsCnt, byrow=TRUE)
			krmm_pckr = raster (krmm_pckr, xmn=region$regionLonMinNominal, xmx=region$regionLonMaxNominal, ymn=region$regionLatMinNominal, ymx=region$regionLatMaxNominal, crs="+proj=lonlat +datum=WGS84")
			writeRaster	(krmm_pckr, paste(DirPrecipKrige, outputDir, "\\", FilePrefixPrecipKrgEstPc, "_", region$regionName, "_", yearKrg, "_", sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE) 
			
			# Output kriged standard error values in percent
			krmm_pcse = matrix (krmm$pcse, region$regionLatCellCntrsCnt, region$regionLonCellCntrsCnt, byrow=TRUE)
			krmm_pcse = raster (krmm_pcse, xmn=region$regionLonMinNominal, xmx=region$regionLonMaxNominal, ymn=region$regionLatMinNominal, ymx=region$regionLatMaxNominal, crs="+proj=lonlat +datum=WGS84")
			writeRaster	(krmm_pcse, paste(DirPrecipKrige, outputDir, "\\", FilePrefixPrecipKrgErrPc, "_", region$regionName, "_", yearKrg, "_", sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
			
		}
		
		if( mode == "t" ) {
			
			# Block kriging of current month, current year, per variogram paramters and prepared data.
			kr = krigblock(region,vgmModel,toLonLat,fromLonLat,prepData$pc[,i],c(CellLonDim,CellLatDim))
			
			# Output kriged prediction values in degrees
			kr_predict = matrix (kr$predict, region$regionLatCellCntrsCnt, region$regionLonCellCntrsCnt, byrow=TRUE)
			kr_predict = raster (kr_predict, xmn=region$regionLonMinNominal, xmx=region$regionLonMaxNominal, ymn=region$regionLatMinNominal, ymx=region$regionLatMaxNominal, crs="+proj=lonlat +datum=WGS84")
			writeRaster	(kr_predict, paste(DirTempKrige, outputDir, "\\", FilePrefixTempKrgEstDg, "_", region$regionName, "_", yearKrg, "_", sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
			
			# Output kriged standard error values in degrees
			kr_krige.var = matrix (kr$krige.var, region$regionLatCellCntrsCnt, region$regionLonCellCntrsCnt, byrow=TRUE)
			kr_krige.var = raster (kr_krige.var, xmn=region$regionLonMinNominal, xmx=region$regionLonMaxNominal, ymn=region$regionLatMinNominal, ymx=region$regionLatMaxNominal, crs="+proj=lonlat +datum=WGS84")
			writeRaster	(kr_krige.var, paste(DirTempKrige, outputDir, "\\", FilePrefixTempKrgErrDg, "_", region$regionName, "_", yearKrg, "_", sprintf("%02d", month), ".bil", sep=""), format="BIL", bandorder="BIL", overwrite=TRUE)
			
		}		
	}
}



prepStnAnmly <- function(mode,region,month) {
	
	# Create station data using 'readclean' function:
	# 'readclean' refers to 'readmon', 'cleanPrecip', and 'cleanTemp'
	# and uses 'region' information.
	stns = readclean(mode,region,month)
	
	# The data has been read per month, selected per years of data,
	# selected per lat/lon in region, and cleaned per std. dev. and NA values
	data = data.frame(stns$data)		# coerce to data frame for ease of using column names
	meta = data.frame(stns$meta)		# coerce to data frame for ease of using column names
	colnames(data) = DataFileColNames
	colnames(meta) = MetaFileColNames
	
	nstn = length(data[,1])
	nyrs = length(data[1,])	
	
	# total number of non-na values per station and per year
	noss = 1:nstn				
	nobs = 1:nyrs   			
	for( i in 1:nstn )  noss[i] = sum(as.numeric(!is.na(data[i,])))
	for( i in 1:nyrs )  nobs[i] = sum(as.numeric(!is.na(data[,i])))	
	
	mn = apply(data,1,mean,na.rm=T)
	sd = apply(data,1,var,na.rm=T)^0.5
	
	# good stations have 25 or more years of observations
	gd = noss >= 25			
	data2 = data[gd,]
	meta2 = meta[gd,]	
	nstn = length(data2[,1])	
	cat(paste("Selected",nstn,"stations having over 25 observations.","\n"))
	
	# plot good stations on map
	par(pty="s") # square plotting region
	plot(meta2$Lon,meta2$Lat,xlab="Longitude",ylab="Latititude",cex=0.5)
	require(maps)
	map(add=T,cex=1.5)
	
	# find column index for beginning and end years for finding anomalies
	# (often use the variogram years, but could be any other appropriate period)
	colYearBegRefMn = which(names(data2)==YearBegRefMn)	
	colYearEndRefMn = which(names(data2)==YearEndRefMn)	
	
	# calculate the mean of each station for finding anomalies
	# from the reference time period (used variogram years, as mentioned above)
	meanRef = apply(data2[,colYearBegRefMn:colYearEndRefMn],1,mean,na.rm=T) 
	
	# calculate anomalies compared to reference mean.
	# prepare matrix for percent anomalies for good stations and all data years
	pc = matrix(NA,nrow=nstn,ncol=nyrs) 
	colnames(pc) = DataFileColNames
	
	# coerce to matrix for correct operation in for loop
	data2 = as.matrix(data2)	
	
	# calculate anomalies for good stations and all data years,
	# in comparison to mean of reference time period
	# (used variogram years, as mentioned above).
	for( i in 1:nstn ) {
		if(mode == "p") {pc[i,] = (data2[i,]+25.)/(meanRef[i]+25.)} 
		if(mode == "t") {pc[i,] = data2[i,] - meanRef[i]}
	}	
	return(list(data=data2,meta=meta2,pc=pc))
}

#
# Clean stations using functions to read and clean.
#
readclean = function(mode,region,month) {
	
	# Create station data using 'readmon' function:
	# 'readmon' and uses 'region' information.
	# The data has been read per month,
	# selected per years in REGION?,
	# selected per lat/lon in region.
	stnsMonYrsReg 	= readmon(mode,region,month)	
	dataMonYrsReg 	= stnsMonYrsReg$data
	metaMonReg 		= stnsMonYrsReg$meta 
	
	# Clean per std. dev. and NA values.
	if( mode == "p" ) {dataMonYrsRegCln  = cleanPrecip(dataMonYrsReg,SigmaPrecip)}
	if( mode == "t" ) {dataMonYrsRegCln  = cleanTemp(dataMonYrsReg,SigmaTemp)}
	
	return(list(data=dataMonYrsRegCln, meta=metaMonReg))
}


#
# Read stations, select by year, and buffer to 10 degrees around region
#
readmon = function(mode,region,month) {  
	
	# Prepare precipitation data/metadata file names
	if( mode == "p" )    {dataFile = paste(DirPrecipData, FilePrefixPrecipData, sprintf("%02d", month), "_", FileSuffixPrecipData, sep="")}
	if( mode == "p" )    {metaFile = paste(DirPrecipMeta, FilePrefixPrecipMeta, sprintf("%02d", month), "_", FileSuffixPrecipMeta, sep="")}
	
	# Prepare temperature data/metadata file names
	if( mode == "t" )    {dataFile = paste(DirTempData, FilePrefixTempData, sprintf("%02d", month), "_", FileSuffixTempData, sep="")}
	if( mode == "t" )    {metaFile = paste(DirTempMeta, FilePrefixTempMeta, sprintf("%02d", month), "_", FileSuffixTempMeta, sep="")}
	
	# Read data and metadata per month	
	dataMon = read.csv(dataFile,header=FALSE)
	metaMon = read.csv(metaFile,header=FALSE)	
	colnames(dataMon) = DataFileColNames
	colnames(metaMon) = MetaFileColNames
	
	# Select years of data
	colYearBegData 	= which(names(dataMon)==YearBegData)
	colYearEndData 	= which(names(dataMon)==YearEndData)
	dataMonYrs 		= dataMon[,colYearBegData:colYearEndData]	
	
	# Select stations within +/- BufferReg (e.g. 10 degrees) of the the region bounding box
	indxStnsInReg 	= 	(metaMon$Lat > (region$regionLatMinNominal-BufferReg)) & 
			(metaMon$Lat < (region$regionLatMaxNominal+BufferReg)) & 
			(metaMon$Lon > (region$regionLonMinNominal-BufferReg)) & 
			(metaMon$Lon < (region$regionLonMaxNominal+BufferReg))	
	dataMonYrsReg 	= dataMonYrs[indxStnsInReg,]
	metaMonReg 		= metaMon[indxStnsInReg,] 
	
	# Count stations before and after selecting per year and region
	numStnsData 	= length(dataMon[,1])
	numStnsReg 		= length(dataMonYrsReg[,1])	
	cat(paste("Read",numStnsData,"stations.","\n"))
	cat(paste("Selected",numStnsReg,"stations within",BufferReg,"degrees of the region.","\n"))
	
	return(list(data=dataMonYrsReg,meta=metaMonReg))	
}

#
# Clean data screens data by a max/min std dev criteria, and fills in missing values
#
cleanPrecip = function(data,sigma) { 
	
	require(MASS)
	cat(paste("Cleaning Data.","\n"))
	
	n   = length(data[,1]) 				# number of selected stations from data
	ny  = length(data[1,]) 				# number of selected years from data
	dobs  = matrix(NA,nrow=n,ncol=ny) 	# matrix of NA values
	
	# Fill in NA matrix with data values
	for( i in 1:n ) {
		dobs[i,] = as.numeric(data[i,])
	}
	
	# Find indices of data that have
	# values below 0 and above 900
	# for assignment of NA values.
	nana  = dobs < 0 | dobs > 900
	
	# Matrix with data and NA values as specified
	dobs[nana]   = NA	
	
	zs = matrix(0,nrow=n,ncol=ny)					# matrix of '0' values, will be z-scores below
	
	colYearBegVgm = which(names(data)==YearBegVgm)	# find column index for beginning year (of variogram?)
	colYearEndVgm = which(names(data)==YearEndVgm)	# find column index for end year (of variogram?)		
	
	for( i in 1:n ) {
		
		d  		= dobs[i,colYearBegVgm:colYearEndVgm]	#select all data in 'variogram' years
		ok 		= !is.na(d)				# indices of cells in range of 'variogram' years that are not NA
		mn 		= mean(d,na.rm=T)		# mean of cells in range of 'variogram' years that are not NA	
		sd 		= var(d,na.rm=T)^0.5	# standard deviation of cells in range of 'variogram' years that are not NA	
		minv 	= min(d,na.rm=T)		# minimum of cells in range of 'variogram' years that are not NA
		maxv 	= max(d,na.rm=T)		# maximum of cells in range of 'variogram' years that are not NA
		
		if( maxv - minv > 600 ) {	# if min to max range of 'variogram' years is > 600 (for this station), 
			zs[i,ok] = sigma+1		# set to zscore outside of sigma parameter, for all years in period of data
		}
		else if( sum(as.numeric(ok)) > 18 && mn > 20 && sd > 20) { #### (e.g. station i=4)
			
			d = d[ok]									# use only non-na (missing) values
			d[d == 0]=1									# set 0 values to 1
			f  = fitdistr(d,"gamma")					## 	Warning ignored by try(): 
			#	'In dgamma(x, shape, scale, log) : NaNs produced'; 
			#	Alleviated by including lower=0, 
			#	but slightly different results for gamma parameters 
			#	'shape', 'scale' and produced a different error on another station row; 
			#	used try() to ignore warning
			
			shape = f$estimate["shape"]
			scale = 1.0/f$estimate["rate"]
			
			d  = dobs[i,]								# select all station data in 'data' years
			ok = !is.na(d)   							# indices of non-na cells in range of 'data' years
			p = pgamma(d[ok],shape,1.0/scale)			# calculate gamma dist probabilities of non-na data in all 'data' years, using gamma dist parameters found from 'variogram' years
			
			p <- ifelse ( p < 0.000001, 0.000001, p) 	# Use 'ifelse' because 'if' only takes a single logical operator, not vectors.
			# (Notice assignment.) p is a vector of station gamma dist probability values of all non-na 'data' years
			
			zs[i,ok] = qnorm(p)							# sets the Normal quantiles of a station using gamma probabilities of all non-na 'data' years
		}
	}
	sigbad = abs(zs) > sigma
	dobs[sigbad] = NA
	cat(paste("Replaced",sum(as.numeric(nana),na.rm=T),"missing values."))
	cat("\n")
	cat(paste("Replaced",sum(as.numeric(sigbad),na.rm=T),"extreme values."))
	cat("\n")
	return(dobs)
}


#
# Clean data screens data by a max/min std dev criteria, and fills in missing values
#
cleanTemp = function(data,sigma) { 
	
	require(MASS)
	cat("Cleaning Data.")
	cat("\n")
	
	n   = length(data[,1])
	ny  = length(data[1,])
	dobs  = matrix(NA,nrow=n,ncol=ny)
	
	for( i in 1:n ) {
		dobs[i,] = as.numeric(data[i,])
	}
	
	nana  = dobs < 0 | dobs > 900
	dobs[nana]   = NA
	
	zs = matrix(0,nrow=n,ncol=ny)
	
	colYearBegVgm = which(names(data)==YearBegVgm)	# find column index for beginning year
	colYearEndVgm = which(names(data)==YearEndVgm)	# find column index for end year		
	
	for( i in 1:n ) {
		
		d  = dobs[i,colYearBegVgm:colYearEndVgm]	
		ok = !is.na(d)
		mn = mean( d,na.rm=T)
		sd = var(  d,na.rm=T)^0.5
		
		if( sum(as.numeric(ok)) > 18 && mn > 20 && sd > 20) {
			d = d[ok]
			d[d == 0]=1
			f  = fitdistr(d,"gamma")
			d  = dobs[i,]
			ok = !is.na(d)         
			zs[i,ok] = (d[ok]-mn)/sd
		}
	}
	sigbad = abs(zs) > sigma
	cat(paste("Replaced",sum(as.numeric(nana),na.rm=T),"missing values.","\n"))
	cat(paste("Replaced",sum(as.numeric(sigbad),na.rm=T),"extreme values.","\n"))
	dobs[sigbad] = NA
	return(dobs)
}



#
# regionLocs creates a 2-d matrix of lat/lon pairs for the background grid.
# The first value is the upper-left corner (i.e. maxlat,minlon).
# The next value is one grid cell to the east.
# The last value will be the lower-right corner (minlat,maxlon).
#
regionLocs = function(region) {
	
	nx = 1 + (region$regionLonMax-region$regionLonMin) / CellLonDim
	ny = 1 + (region$regionLatMax-region$regionLatMin) / CellLatDim
	nxny = nx * ny
	cat(paste("Number of longitudinal values in region: ", nx,"\n"))
	cat(paste("Number of latitudinal values in region: ", ny,"\n"))
	cat(paste("Number of kriging output locations in region: ", nxny,"\n"))
	
	locs = matrix(nrow=nxny,ncol=2)
	for( y in 1:ny ) {
		for( x in 1:nx ) {
			idx = ((y-1)*nx) + x 
			locs[idx,1] = region$regionLonCellCntrs[x]
			locs[idx,2] = region$regionLatCellCntrs[y]
		}
	}
	colnames(locs) = c("lon","lat")	
	return(locs)
}


#
# Read in 0.25 degree monthly 'RefData' climatology 
# (i.e. FCLIM) as reference data to find anomalies.
#
readRefData <- function(month) {
	
	require(raster)
	
	#Specify path for data
	inputFile = paste(DirRefData, FilePrefixPrecipRefData, sprintf("%02d",month), FileSuffixPrecipRefData, sep="")
	
	#read the data
	rasterData <- raster(inputFile) 
	return(rasterData)
	print("Extent of Reference Data")
	print(extent(rasterData))
	
}


###################################
# Define function 'krigblock' to 
# perform block kriging interpolation
# using region, variogram parameters,
# locations to krige from and to,
# the prepared data and the
# block (raster cell) dimensions.

krigblock <- function(region,vgmModel,toLocs,fromLocs,data,blocks) {
# Example usage below; note data is square root of percent anomalies; block definition in Lon/Lat dimensions
# kr = krigblock(region,vgmModel,toLonLat,fromLonLat,pc[,i]^0.5,c(CellLonDim,CellLatDim)) 
	
	# Load R package
	require(gstat)
	
	# Prepare variables for kriging ...
	obs  	<- as.numeric(data)
	n    	<- length(obs)
	
	tmp 	<- matrix(NA,nrow=n,ncol=3)
	tmp[,1] <- fromLocs[,1] 
	tmp[,2] <- fromLocs[,2]
	tmp[,3] <- obs
	
	# 'De-duplicate': remove duplicates, replace their values with mean of duplicates	
	tmp <- dedup(tmp)						
	
	# Identify stations without NA values in data (square root of percent anomalies)
	gid <- !is.na(tmp[,3])					
	tmp <- tmp[gid,]						
	
	# Set header labels for prepared data and kriged locations
	colnames(tmp)  		<- c("lon","lat","obs")
	colnames(toLocs) 	<- c("lon","lat")
	
	# Coerce data types
	tmp    <- data.frame(as.matrix(tmp))
	toLocs <- data.frame(as.matrix(toLocs))
	
	# Alternative kriging function below for other purposes.
	# t <- krige(obs~lon+lat,~lon+lat,data=tmp,model=vgmModel,newdata=tolocs,nmax=20,block=blocks)
	t <- krige(obs~1,~lon+lat,data=tmp,model=vgmModel,newdata=toLocs,nmax=KrigeNMax,block=blocks)
	
	# Create list for returning kriging predictions and standard errors
	kr <- list(predict=t[,3],krige.var=t[,4])
	
	return(kr)
}




dedup = function(data) { 							# "deduplicate": remove duplicate lon, lat values
	
	lon  = data[,1]
	lat  = data[,2]
	ll   = cbind(lon,lat)
	lln = unique(ll)								# return only unique lon, lat data
	dn  = matrix(NA,nrow=nrow(lln),ncol=ncol(data))	# prepare matrix for output
	colnames(dn) = colnames(data)
	nc = ncol(data)
	for( i in 1:nrow(dn) ) {
		sub = lon == lln[i,1] & lat == lln[i,2]		#identify original rows (stations) with the same lon, lat values as the unique lon, lat
		for( j in 1:nc )
			dn[i,j] = mean(data[sub,j])				#for each column, place the mean of each row or duplicate sets of rows into ouput data
	}
	return(dn)
}





#
# pct2mm rescales pct estimates and percent errors by multiplying them by their 
# background reference data (i.e. FCLIM, long term average) values.
#
pct2mm <- function(region,ltm,pr.pc,pr.er,locs) {
# Example usage below
# pct2mm(region,RefData,pr.pc,pr.er,toLonLat)	
	sclX  <- 1.0 / CellLonDim	
	sclY  <- 1.0 / CellLatDim	
	lon  <- locs[,1]
	lat  <- locs[,2]
	
	ullat <- locs[1,2]
	ullon <- locs[1,1]
	
	#### Half-cell shift for addressing 
	#### the raster() reading the extents of the file
	#### compared to the kriged cell center,
	####
	#### Plus  half-cell adjust. (e.g.  0.125 ullon + 0.25 Cell/2 =  1 MinLonBG)
	#### plus works for negatives: if ullon= -179.875, then x1=   1 and 'row' cols = 1:240
	#### plus works for positives: if ullon=    0.125, then x1= 721 and 'row' cols = 721:960
	x1   <- ((ullon + CellLonDim/2) - MinLonBG) * sclX 
	
	#### Plus  half-cell adjust. (e.g. 49.875 ullat + 0.25 Cell/2 = 50 MaxLatBG) ### For "Region_0_50N_0_60E", y1=   0; this seems to make sense for the 'row'
	#### plus works for positives: if ullat=   49.875, then y1=   0 and 'row' row =   1
	#### plus works for negatives: if ullat=   -0.125, then y1= 200 and 'row' row = 201
	y1   <- (MaxLatBG - (ullat + CellLonDim/2)) * sclY 	
	
	mmkr <- matrix(NA,nrow=nrow(locs))
	mmse <- matrix(NA,nrow=nrow(locs))
	pckr <- matrix(NA,nrow=nrow(locs))
	pcse <- matrix(NA,nrow=nrow(locs))
	pmm  <- matrix(NA,nrow=region$regionLonCellCntrsCnt)
	emm  <- matrix(NA,nrow=region$regionLonCellCntrsCnt)
	
	n <- 1
	for( y in 1:region$regionLatCellCntrsCnt ) {
		
		row <- ltm[(y1+y),x1:(x1+region$regionLonCellCntrsCnt-1)]
		p   <- pr.pc[n:(n+region$regionLonCellCntrsCnt-1)]
		e   <- pr.er[n:(n+region$regionLonCellCntrsCnt-1)]^0.5
		
		id  <- row >= MinValBG & row < MaxValBG & !is.na(p) & !is.na(row)
		pmm[id]  <- row[id] * p[id]
		emm[id]  <- row[id] * e[id]
		pmm[!id] <- NA
		emm[!id] <- NA
		p[!id] <- NA
		e[!id] <- NA
		
		mmkr[n:(n+region$regionLonCellCntrsCnt-1)] <- pmm
		mmse[n:(n+region$regionLonCellCntrsCnt-1)] <- emm
		pckr[n:(n+region$regionLonCellCntrsCnt-1)] <- p
		pcse[n:(n+region$regionLonCellCntrsCnt-1)] <- e
		
		n <- n + region$regionLonCellCntrsCnt
	}
	krmm <- list(mmkr=mmkr,mmse=mmse,pckr=pckr,pcse=pcse)
	
	return(krmm)
}



#
# run_krgMosaicPrecip
# for running 'krgMosaic' on a multiple precipitation image sets
# per directories of months, and files of
# years, output types, and regions
#
run_krgMosaicPrecip <- function (){
	krgMosaic 		(mode="p", DirTagMonthsKrg, FileTagYearsKrg, FileTagKrgOutputPrecip, FileTagRegions)
}

#
# run_krgMosaicTemp
# for running 'krgMosaic' on a multiple temperature image sets
# per directories of months, and files of
# years, output types, and regions
#
run_krgMosaicTemp 	<- function (){
	krgMosaic 		(mode="t", DirTagMonthsKrg, FileTagYearsKrg, FileTagKrgOutputTemp, FileTagRegions)
}


#
# krgMosaic
# for mosaicing multiple image sets
# per directories of months, and files of
# years, output types, and regions
#
krgMosaic <- function(mode, dirTagMonthsKrg, fileTagYearsKrg, fileTagKrgOutput, fileTagRegions) {
	
	require(raster)
	
	for (mon in 1:length(dirTagMonthsKrg)) {
		
		if( mode == "p" ) {
			outputDir = paste(sprintf("%02d", dirTagMonthsKrg[mon]),"\\",sep="")		
		}
		
		if( mode == "t" ) {
			outputDir = paste(sprintf("%02d", dirTagMonthsKrg[mon]),"\\",sep="")		
		}
		
		for (krgYear in 1:length(fileTagYearsKrg)){
			
			for (krgOut in 1:length(fileTagKrgOutput)){
				
				inputFile = paste(DirPrecipKrige, outputDir, fileTagKrgOutput[krgOut], "_", fileTagRegions[1], "_", fileTagYearsKrg[krgYear], "_", sprintf("%02d", dirTagMonthsKrg[mon]), ".bil", sep="")
				mosaicFile = raster(inputFile)
				
				for (reg in 2:length(fileTagRegions)){
					inputFile 	= paste(DirPrecipKrige, outputDir, fileTagKrgOutput[krgOut], "_", fileTagRegions[reg], "_", fileTagYearsKrg[krgYear], "_", sprintf("%02d", dirTagMonthsKrg[mon]), ".bil", sep="")
					region 		= raster(inputFile)
					mosaicFile 	= mosaic(mosaicFile,region,fun=mean) # fun=mean averages overlapping values
				}
				
				outMosaicName 	= paste(DirPrecipKrige, outputDir, fileTagKrgOutput[krgOut], "_",fileTagYearsKrg[krgYear],"_",sprintf("%02d", dirTagMonthsKrg[mon]),".bil",sep="")
				projection	(mosaicFile) <- "+proj=lonlat +datum=WGS84"
				writeRaster	(mosaicFile, filename=outMosaicName, format="BIL", bandorder="BIL", overwrite=TRUE)	#Support .hdr and .stx files automatically written.
			}
		}
	}
}




#
# run_krgLandMaskPrecip
# for running 'krgLandMask' on a multiple precipitation image sets
# per directories of months, and files of
# years, and output types
#
run_krgLandMaskPrecip 	<- function(){
	krgLandMask 		(mode="p", DirTagMonthsKrg, FileTagYearsKrg, FileTagKrgOutputPrecip)
} 

#
# run_krgLandMaskTemp
# for running 'krgLandMask' on a multiple temperature image sets
# per directories of months, and files of
# years, and output types
#
run_krgLandMaskTemp 	<- function(){
	krgLandMask 		(mode="t", DirTagMonthsKrg, FileTagYearsKrg, FileTagKrgOutputTemp)
} 

#
# krgLandMask
# for masking multiple image sets
# per directories of months, and files of
# years, and output types
#
krgLandMask <- function (mode, dirTagMonthsKrg, fileTagYearsKrg, fileTagKrgOutput) {
	
	require(raster)
	
	for (mon in 1:length(dirTagMonthsKrg)) {
		
		if( mode == "p" ) {
			outputDir = paste(sprintf("%02d", dirTagMonthsKrg[mon]),"\\",sep="")	
		}
		
		if( mode == "t" ) {
			outputDir = paste(sprintf("%02d", dirTagMonthsKrg[mon]),"\\",sep="")		
		}
		
		for (krgYear in 1:length(fileTagYearsKrg)){
			
			for (krgOut in 1:length(fileTagKrgOutput)){	
				
				inputImageName = paste(DirPrecipKrige,outputDir,fileTagKrgOutput[krgOut],"_",fileTagYearsKrg[krgYear],"_",sprintf("%02d", dirTagMonthsKrg[mon]),".bil",sep="")
				inputImage = raster(inputImageName)
				
				maskImageName = paste(DirLandMask, FileLandMask, sep="")
				maskImage = raster(maskImageName)
				maskImage[maskImage == 0] <- NA
				maskImage = crop(maskImage, inputImage) # Crop maskImage to the extent of the inputImage
				
				outputImageName = paste(DirPrecipKrige,outputDir,fileTagKrgOutput[krgOut],"_",fileTagYearsKrg[krgYear],"_",sprintf("%02d", dirTagMonthsKrg[mon]),"_land",".bil",sep="")
				outputImage = mask(inputImage, maskImage)
				projection	(outputImage) <- "+proj=lonlat +datum=WGS84"
				writeRaster (outputImage, filename=outputImageName, format="BIL", bandorder="BIL", overwrite=TRUE)
			}
		}
	}
}


